/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TABound
  \brief To calculate bound state radial wavefunction of a valence nucleon in
  central potential (Vc+VN+VLS) of a nucleus.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2020/07/08
  \date Last modified: 2020/09/20 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#ifndef _TABound_h_
#define _TABound_h_

#include <string>
#include <vector>
#include "TAYaml.h"
#include "TAODESolver.h"

using std::string;
using std::vector;

class TABound : public TAODESolver{
public:
  /// \param: v: input parameters
  TABound(const TAYaml &v);
  virtual ~TABound();

  double V(double r) const; ///< \retval V(r) = (V0+VSO+VC)*f2MuOverHbar2
  // potential for the unit-test method BoundExample
  double VT(double r) const; ///< \return total potential: V(r)+VL (VL is the centrifugal barrier)
  int Getl() const{ return l; }
  /// \retval Note that Rl=u*r; u = Rl/r. This method returns Rl, i.e. u*r, NOT u
  double u(double r); ///< u(r)=r*Rl
  double Rl(double r);
  double E() const{ return fE; }

  /// the ODE set: du/dr=u1; du1/dr=-(lambda+v)*u; du2/dr=0,  where E=hbar^2/(2mu)*lambda ///
  /// the two point border condition: u(0)=0; du(0)/dr=1; u(\infty)=0; ///
  /// the ODE set itself: dyi/dx = f_i(x,y1,..yN)
  virtual void derivs(double x, const double *y, double *dydx) override;
  /// calculates the n-vector y[0..n-1] (satisfying the starting boundary conditions, of course)
  /// given the freely specifiable variables of v[0..n2-1] at the initial point x1
  // solve the radial wavefunction //
  void Bound(); ///< for solving the objective radial equation
  void BoundIterE(); ///< iterate Bound by adjusting fV0 to reproduce S_N*
  void BoundIterR(); ///< iterate Bound by adjusting fR0,fRC,fRS to reproduce r_rms
  /// iterate BoundIterE by adjusting fV0 to reproduce S_N*, and then adjusting fR0, to reproduce r_rms
  void BoundIterER();
  void BoundIterNewton(); ///< BoundIterER with Newton global root-finding approach
  void BoundIter(); ///< self-adaptive BoundIter according to availabilities of fRref and fRref
  bool HasSolved(){ return fSolved; } ///< if Bound() has been called or not
  const char *Orbit(){ return fOrbit.c_str(); }

  static double ZERO; // fx1
  static double INFTY; // fx2
private:
  bool fSolved; ///< if Bound() has been called or not

  double fZc, fAc; ///< identity of the core
  double fZv, fAv; ///< identity of the valence nucleon
  /// single-particle state for the valence nucleon
  int n, l;
  double j;
  /// the nuclear potential: WS form
  double fV0, fR0, fA0; ///< central NN force
  double fVS, fRS, fAS; ///< central NN force
  double fRC; ///< charge radius for Coulomb potential

  double fCoeVS, fCoeCoul, fCoel; ///<  VS*RPICOMP^2/AS*<s.l>, Zc*Zv*e2, hbar^2/(2*mu)*l(l+1)
  double f2MuOverHbar2; // 2mu/hbar^2

  double *fxx, *fyy2; ///< the concatenated x and y^2
  int fcnt; // the length of the solution fyy
  double fE; ///< the solved eigen energy in MeV
  double fRrms; ///< the solved rms radius in fm

  double fEref, fRref; ///< reference separation energy and HF r_rms
  int fFindOpt; ///< 0: Newton global search; 1: RE search: E embedded in R
  string fOrbit; ///< the orbit info
};

#endif
